Nature Conservancy units

library(dplyr)
library(here)
library(mapview)
library(tigris)
library(sf)

mapviewOptions(basemaps.color.shuffle = FALSE)

ca_blocks_clean <- readRDS(here("Analysis", "Data", "Blocks", "v202506", "CA_blocks_v202506.RDS")) %>%
  filter(WATER_SLIVER == 0 & CA_SLIVER == 0)
cntr_crds <- c(mean(st_coordinates(ca_blocks_clean)[, 1]),
               mean(st_coordinates(ca_blocks_clean)[, 2]))

# Load State Park data
tnc <- st_read(here("Analysis", "Data", "CPAD", "cpad_2024a", "CPAD_2024a_Holdings.shp")) %>%
  filter(AGNCY_NAME == "The Nature Conservancy" | MNG_AGNCY == "The Nature Conservancy") %>%
  st_transform(st_crs(ca_blocks_clean)) %>%
  st_make_valid() 
  #%>%
  #st_cast("POLYGON") %>%
  #mutate("temp_ID" = row_number(.))

tnc_temp <- tnc %>%
  mutate(area = st_area(.) %>% units::drop_units()) %>%
  st_drop_geometry() %>%
  group_by(UNIT_NAME) %>%
  summarize(unitname_area =sum(area))

# Intersect with CA-BBA blocks, count number of parks by block, and record thier IDs and the ID of park with largest footprint in block
xFeatures <- ca_blocks_clean %>% select(BLOCKNAME)
yFeatures <- tnc %>% select(HOLDING_ID, ACCESS_TYP, UNIT_ID, UNIT_NAME, AGNCY_NAME, MNG_AGNCY)
intersections <- st_intersects(x = xFeatures, y = yFeatures)

library(progress)
pb <- progress_bar$new(format = "[:bar] :current/:total (:percent)", total = dim(xFeatures)[1])

library(purrr)
tnc_intersect <- map_dfr(1:dim(xFeatures)[1], function(ix){
  pb$tick()
  st_intersection(x = xFeatures[ix,], y = yFeatures[intersections[[ix]],])
})

tnc_intersect_details <- tnc_intersect %>%
  mutate(area = st_area(geometry)) %>%
  as.data.frame() %>%
  st_drop_geometry() %>%
  select(BLOCKNAME, UNIT_NAME, area) %>%
  mutate(area = units::drop_units(area)) %>%
  group_by(BLOCKNAME, UNIT_NAME) %>%
  left_join(tnc_temp %>% st_drop_geometry()) %>%
  summarize(unit_area = sum(area),
            unit_pct_area_in_block = (sum(area) / mean(unitname_area) * 100) %>% round(5),
            .groups = 'drop')
tnc_intersect <- tnc_intersect_details %>%
  group_by(BLOCKNAME) %>%
  summarize(
   TNC_UNIT_COUNT = length(unique(UNIT_NAME)),
   TNC_UNIT_PRIMARY = UNIT_NAME[which.max(unit_area)],
   TNC_UNIT_MOST_LOCALIZED = UNIT_NAME[which.max(unit_pct_area_in_block)],
   TNC_UNIT_MOST_LOCALIZED_PCT = unit_pct_area_in_block[which.max(unit_pct_area_in_block)],

   TNC_UNIT_ALL = paste(unique(UNIT_NAME), collapse = "; "),
    .groups = 'drop'
  )

# Join State Park info back to full blocks
ca_blocks_clean <- ca_blocks_clean %>%
  left_join(tnc_intersect, by = "BLOCKNAME") 

# Calculate number of blocks with TNC holdings
n_blocks <- ca_blocks_clean %>% nrow()
n_blocks_TNC <- ca_blocks_clean %>% filter(!is.na(TNC_UNIT_MOST_LOCALIZED)) %>% nrow()